********************************************************************************
***********************Columns (1) and (2)**************************************
********************************************************************************

clear all
set more off
global dpath "[Userpath]\data and code"

use "$dpath/data/countyEmission.dta", clear

tostring countycode,g( countycode_str )
g citycode_str=substr( countycode_str ,1,4) 
g a_str=substr( countycode_str ,-2,. )
destring a_str,g(a) 

gen exempt = 1 if citycode_str=="1521" | 									    ///
					  citycode_str=="2224" | 								    ///
					  citycode_str=="4190" | 								    ///
					  citycode_str=="4521" | 								    ///
					  citycode_str=="4600" | 								    ///
					  citycode_str=="5323" |								    ///
					  citycode_str=="5325" | 							        ///
					  citycode_str=="5328" | 								    ///
					  citycode_str=="5329" | 								    ///
					  citycode_str=="5331" | 								    ///
					  citycode_str=="6541" | 								    ///
					  citycode_str=="6540" |								    ///
					  citycode_str=="6543"
replace a_str="00" if a<20 & exempt != 1
g countycode2_str= citycode_str+ a_str
destring countycode2_str,g( countycode2 )
drop countycode_str a_str citycode_str a countycode2_str exempt
order countycode2,after( countycode )


sort countycode2 cic2 year
foreach var in wastewater cod nh3n gas so2 dust {
	bysort countycode2 cic2 year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
}
duplicates drop countycode2 cic2 year, force

rename countycode countycode_ini
rename countycode2 countycode
drop dq
order countycode countycode_ini


***************************** match with station info **************************
joinby countycode using "$dpath/data/countyStation100km.dta", unmatched(both)
tab _merge

keep if _merge == 3
drop _merge

gen citycode = int(countycode/100)
gen provincecode = int(citycode/100)

joinby countycode year using "$dpath/data/countyControls100km.dta", unmatched(both)
tab _merge

keep if _merge == 3
drop _merge


joinby countycode using "$dpath/data/auto_up.dta", unmatched(both)
tab _merge

drop if _merge == 2
drop _merge
replace autoyear_up = 33000 if autoyear_up == .

gen auto_up = 1 if year>=autoyear_up
replace auto_up = 0 if year<autoyear_up

gen llgdpC = log(lgdpC)
drop if llgdpC == .

gen llgdp = log(lgdp)
drop if llgdp == .

foreach var in wastewater cod nh3n gas so2 dust {
	bysort countycode year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
	gen l`var' = log(`var')
}
duplicates drop countycode year, force

xtset countycode year

winsor2 lwastewater, replace cuts(0.5 99.5) trim
winsor2 lcod, replace cuts(0.5 99.5) trim
winsor2 lgas, replace cuts(0.5 99.5) trim
winsor2 lso2, replace cuts(0.5 99.5) trim
winsor2 ldust, replace cuts(0.5 99.5) trim

xtset countycode year

gen manualupstream = 1 if manualyear != 33000
replace manualupstream = 0 if manualupstream == .

gen treated = 1 if autoupstream == 1 & year >= autoyear
replace treated = 0 if treated == .

gen treatedM = 1 if manualupstream == 1 & year >= manualyear
replace treatedM = 0 if treatedM == .

cd "$dpath\results"

reghdfe lwastewater treated, absorb(countycode citycode#year) cluster(citycode)
estimates store reg1
reghdfe lwastewater treated treatedM llgdpC, absorb(countycode citycode#year) cluster(citycode)
estimates store reg2
esttab reg1 reg2 using TableAVIII.rtf,replace b(%9.3f) se scalar(N r2_a) title(Columns 1 and 2)  star(* 0.1 **  0.05  *** 0.01)	 keep(treated)
est clear

********************************************************************************
***********************Columns (3) to (6)***************************************
********************************************************************************
set     more off
set     matsize 800
clear   all
global  dpath = "G:\水污染与健康\data and code"

use "$dpath\data\deathratecontrol.dta",clear

cd "$dpath\results"

reghdfe totaldigestivedr auto100, absorb(code year)vce (cluster citycode)
estimates store reg1
reghdfe totaldigestivedr auto100 manual100 lngdp lnbed, absorb(code year)vce (cluster citycode)
estimates store reg2
winsor2 lifeexpectancy, replace cuts(1 99) trim
reghdfe lifeexpectancy auto100, absorb(code year)vce (cluster citycode)
estimates store reg3
reghdfe lifeexpectancy auto100 manual100 lngdp lnbed, absorb(code year)vce (cluster citycode)
estimates store reg4
esttab reg1 reg2 reg3 reg4 using TableAVIII.rtf,append b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Columns 3 to 6) keep(auto100)
est clear